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Abstract. A set of bi-orthogonal potential-density basis functions is introduced to model the 
density and its associated gravitational field of three dimensional stellar systems. Radial compo- 
^— «j . nents of our basis functions are weighted integral forms of spherical Bessel functions. We discuss 
' on the properties of our basis functions and demonstrate their shapes for the latitudinal Fourier 
, number I — 2. 
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CN " 1. Introduction 

,—1 ^ Potential-density basis sets are used for expanding the solution of Poisson's equation 
. during N-body simulations and stability analysis of stellar systems. Spherical harmonics 
are the most efficient functions for expansions in the azimuthal and latitudinal directions 
Q . but finding bi-orthogonal components in the radial direction, which are expressible in 
^ ' terms of ordinary or special functions, is not a straightforward task. A numerical method 
^ . has been introduced by Weinberg (1999) for constructing the basis functions in the radial 
I I ■ direction. His method always gives a bi-orthogonal set whose components oscillate in the 
radial direction. However, the amplitude of oscillations and the locations of zeros of such 
' functions do not necessarily assure a good performance. Especially, when one of the peaks 
, of radial components is higher than others [this happens for the Clutton-Brock (1973) 
— . ] and Hernquist & Ostriker (1992) functions], computations are rapidly biased. Here we use 
Glutton-Brock's (1972) idea and construct a basis set using the weighted integral form 
of spherical Bessel functions. The functions that we find have shown a good performance 
, in the stability analysis of spherical systems ( Jalali & Hunter 2007) . 
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2. A new potential density pair 

■ The Laplace operator is Hermitian, and therefore, its eigenfunctions always form a 
complete and bi-orthogonal set (Arfken 1985). A three dimensional generalization of 
Glutton-Brock's (1972) method suggests to introduce 
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as the potential and density basis functions. Here gni{k) are some orthogonal functions 
and a is a length scale. ji{r) and YimiO, (j)) are the spherical Bessel functions and spherical 
harmonics, respectively. Their combination ji(r)y;m(0, 0) satisfies Poisson's equation. It 
is straightforward to show that V^$/mri = 4:TrGpimn- Denoting as Laugerre functions 
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and considering their orthogonahty, the choice of gin{k) 
following bi-orthogonality relation between ^imn and pimn 



¥L^^{2k)e^'' leads to the 
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The above potential and density functions oscillate in the r-dircction and the amplitude 
of oscillations decays slowly. Peaks of these basis functions accumulate near the center 
by increasing n, and therefore, they are suitable for stability analyses that utilize matrix 
methods. The first potential fimction in the radial direction is ^Qoir) = — r~^arctanr 
whose profile falls off similar to as r ^ oo. Our basis functions can be computed 
both numerically and analytically. Since the spherical Bessel functions and their combina- 
tions with Laugerre fimctions decay rapidly, definite integrals in (2.1) and (2.2) converge 
very fast by adopting standard numerical schemes. For instance, we have computed and 
displayed in Figure 1 the radial potential and density basis functions for / = 2 and sev- 
eral values of n. We have also examined the performance of our basis set in modeling 
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Figure 1. Potential (left) and density (right) basis functions for I = 2. Solid, dashed, dash-dotted 
and dotted lines correspond to n =0, 1, 2, and 3, respectively. Densities have been scaled to 
their maximum values. 

the spherical and oblate Kuzmin-Kutuzov models. Our experiments show that taking a 
handful of expansion terms reconstructs the original model with a relative error less than 
%1. 

3. Conclusions 

Our basis functions constitute a bi-orthogonal set over the domain ^ r < oo. By 
adjusting the length scale a, it is possible to model stellar systems of different core radii 
and spatial extensions. The wave length of our radial basis functions decreases towards 
the center while their wave amplitudes are not changed substantially. These properties 
are useful for the instability analysis of three dimensional stellar systems, for instabilities 
occur in central regions of stellar systems and accurate calculation of perturbed physical 
quantities demands high-resolution functions near the center. 
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